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Abstract 

Statistical methods for the extraction of a small shift in broad data distributions are examined by means of 
Monte Carlo simulations. This work was originally motivated by the CERN neutrino beam to Gran Sasso (CNGS) 
experiment for which the OPERA detector collaboration reported in a broad distribution a time shift by ±7.1 ns 
despite the fact that the fluctuation of the average time is with ±23.8 ns much larger. Although the physical result 
has been withdrawn, statistical methods that make such an identification in a broad distribution possible remain 
of interest. 
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1. Introduction 

In highly publicized CERN announcements [1] it 
was claimed that neutrinos from the CNGS arrived 
at Gran Sasso 



St 



-57.8 ±7.8 (stat. 



-8.3 
-5.9 



(sys.) 



(1) 



too early, violating the St = limit set by the 
speed of light. Meanwhile, it appears that initially 
overlooked systematic errors have wiped out the 
result, but the identification of the statistical part 
remains of interest as it exemplifies the extraction 
of a small shift from a broad distribution. The pur- 
pose of this article to shed light on subtleties of an 
analysis, which leads to the statistical part of the 
estimate (1). 

The CNGS sample of 15 223 neutrinos was pro- 
duced in extractions that last about 10,500 ns each. 



Two different types of extractions were used lead- 
ing to probability densities (PD) 



p k (t), fc = l,2 



(2) 



for neutrinos departure times, which are repro- 
duced here in Fig. 1 within a precision allowed by 
reading them off from Fig. 11 of [1]. The PD used 
in our paper have been discretized in intervals of 
1 ns and can be downloaded from the author's web- 
site [2]. 

One can now perform a statistical bootstrap [3] 
analysis by Monte Carlo (MC) generation of de- 
parture times with the PD of Fig. 1. This is already 
remarked in [1] , where the application remains lim- 
ited to testing of their maximum likelihood proce- 
dure on a sample of 100 MC data sets. As the MC 
generation of departure times can be repeated al- 
most arbitrarily often with distinct random num- 
bers, one can analyze and verify statistical mcth- 
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Fig. 1. Departure time probability densities modeled after Fig. 11 of Ref. [1]. 



ods that one wants to apply to discover a shift in 
the data. 

For the uniform PD over 10 500 ns, it has been 
noted [4] that with n = 16 111 events the variance 
of the departure time average t is approximately 
At = 24ns, i.e., much larger than the statistical 
error bar in Eq. (1). It will be discussed in this pa- 
per that the time shift St (1) defined in [1] behaves 
indeed quite differently than a statistical fluctua- 
tion of the time average 



St 



rti Sti + ri2 St2 



St: 



M3) 



Here U are the measured departure time averages, 
tj are the mean departure times obtained from 
the corresponding PD, and rii are the numbers of 
events in each extraction. The distinction between 
St (1) and St (3) is made by an overline on t or not. 
Obviously, 



(St) = (St) 



(4) 



holds for the expectation values, but their error 
bars behave differently. 

To set the groundwork, it is shown in section 2 for 
the uniform distribution that a shift St = —57.8 ns 
can be identified with certainty (probability to miss 
it < 10 -36 ) when there are 15 223 events and the 
departure time range is 10 500 ns. In section 3 the 
MC generation of departure times is described. 
Section 4 gives examples of descriptive histograms 
from MC data. Suggested by the uniform distribu- 



tion, the front tails of the distributions are of par- 
ticular interest. For their study the cumulative dis- 
tribution function (CDF) is far better suited than 
histograms, because it allows easily to focus on out- 
liers. This is investigated in section 5. To estimate 
the shift value St, the maximum likelihood method 
is used in [1]. In section 6 features of this method 
are calculated by applying it to a large number of 
MC generated departure time samples. 

Independently of the special example, the ap- 
proaches discussed in sections 2 to 6 are of interest, 
because they address the general problem of ex- 
tracting a precise estimate of a shift from a broad 
distribution. Summary and conclusions follow in 
section 7. 



2. Uniform Distribution 

Using the uniform PD over a time window of 
10 500 ns, the standard deviation of the average 



1 n 
n t-f 



(5) 



is for n = 15 223 events much larger than the sta- 
tistical error bar quoted in Eq. (1), namely approx- 
imately 



At = 25 ns. 



(6) 



2 



1.2 
1 

0.8 
0.6 
0.4 
0.2 




Events 
Impossible 

-5000 



5000 
t[ns] 



Missing 
Events 

10000 15000 




0.7424 



2010 2012 2014 2016 2018 2020 2022 2024 
t[ns] 



Fig. 2. Uniform distribution: Impossible (left) and missing 
events (in the enlarged thickness of the right border) . 

How can this be? That the average (5) fluctuates 
with the variance (6) is unavoidable. However, the 
effect we are after is a systematic shift of each de- 
parture time by an amount St = —57.8 ns. Again 
for the uniform uniform distribution, drawn in 
Fig. 2, it is easily illustrated that this can very 
well be identified. Events indicated on the left of 
the figure are impossible unless there is a shift. 
Now, with a shift of —57.8 ns the probability to 
find a particular event to the left of the uniform 
PD is given by 

p = 57.8/10 500 = 0.005505... (7) 

and the probability to find none is 

(1-p) 15672 = HT 36 - 5 . (8) 

The distance of the smallest time from the left edge 
of the uniform PD is a lower bound on St and a 
direct estimate for the time shift (1) is 

St = ni eft 10500ns/15223, (9) 

where ni e ft is the number of events observed on the 
left outside of the uniform PD. Confidence limits 
can be established from the binomial distribution. 

For the tiny range of 57.8 ns indicated by the 
somewhat thicker line on the right side of the uni- 
form PD in Fig. 2, the situation is the other way 
round. It has to be empty when there is a shift by 
St = —57.8 ns. The probability that this happens 
by chance when there is in fact no shift is also given 
by (9). When St is not known the distance of the 



Fig. 3. Enlargement of the approximated PD for model 1 
over a small time region. 

largest measured time from the right edge of the 
uniform PD is an upper bound on St. 

We do not pursue the uniform PD any further, 
because we are interested in the more complicated 
case of the less sharp PD of Fig. 1. 



3. MC generation of departure times 

As mentioned in the introduction, the PD of 
Fig. 1 have been discretized in 1 ns intervals and 
are available on the Web [2] . The resolution of 1 ns 
allows for easy MC generation of departure times 
and is sufficient for the intended accuracy of the es- 
timate of a shift. In the following our thus defined 
models are labeled by k — 1,2. The probabilities 
as function of time t are defined by 

Pk{t) = p fc (**) for i* < t < i l + 1 , (10) 

where i are integer times in ns units. As it is con- 
venient for the MC generation of departure times, 
the normalization for the discretized PD is (dis- 
tinct from Fig. 1) chosen so that 

P r x = max[p fc (i)] = 1 (11) 

holds. Proper normalizations X^PfcW = 1 
could still be achieved by choosing instead of ns 
some unconventional unit for At^. For the gener- 
ation of correctly distributed random times this is 
irrelevant. Model 1 probabilities p\ (t) are enlarged 
in Fig. 3 for a short time range. 
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Fig. 4. Typical histograms of departure times together with the expectation from the PD of Fig. 1. For the (online) green 
entries the generated time distribution is shifted by 57.8 ns to the left. 

After discretization the smallest z™ n and largest i\ = (t)i = 5942.2ns, t 2 = (t)2 = 6658.7ns ,(12) 



times with non-zero values are 



if" 1 = 359 and 



= 12 and 



= 15 368, 
= 19 877. 



These ranges include a number of zero probabili- 
ties in particular for the large i l values. MC gen- 
erated departure times t 3 k (J = 1, . . . , n) with n 
the number of data are obtained from uniformly 
distributed random numbers t in a range enclos- 
ing (ir 11 ,^ + !). here chosen to be (1,22000): 
for model k a proposed random time t is accepted 
with probability for «* < t < i* + l. If t is re- 

jected, the procedure is repeated until a value gets 
accepted, which then becomes a data point t 3 k . Ac- 
ceptance rates were close to 40% and the random 
number generator of Ref. [5] has been used. To 
avoid direct hits of i™ m and associated rounding 
problems, the proposed random times are shifted 
by +2~ 25 , which is half the discretization [6] of 
these random numbers. 



4. Histograms 

In this section we pursue descriptive statistics 
and show histograms from samples of n = 7 612 
MC generated departure times per model. The ex- 
act departure time mean values for our discretized 
PD of Fig. 1 are (the hat indicates exact): 



where model 1 and 2 are again labeled by the corre- 
sponding subscripts. We also know from these PD 
the exact standard deviations for time averages ti 
and t 2 , each over n = 7 612 events: 

Ati = 33.4 ns, Ai 2 = 33.8 ns. (13) 
Their combined standard deviation is 



At = 23.8 ns. 



(14) 



which agrees almost with the estimate (6) from the 
uniform PD. 

Together with the PD, Fig. 4 shows (red online) 
histograms of 150 ns bin width for departure times 
from a typical [7] MC generated sample of 7612 
events per model and again the same histograms 
(green online) shifted by —57.8 ns. This illustrates 
how small the effect is, which ought to be unam- 
biguously identified. 

Rounded to nearest integers the mean values for 
the time averages of the data of Fig. 4 are 



h = (5936 ± 34) , ns, (model 1) 
t 2 = (6703 ± 34) , ns, (model 2) 



(15) 
(16) 



Both are within their error bars consistent with the 
mean values (12) of the PD. The actual statistical 
fluctuations of these two MC data sets are (using 
h = 5935.5 and t 2 = 6703.2) 



4 





Fig. 6. Model 1: Zoom of the leading (left plots) and trailing (right plots) MC generated departure time histograms (red 
online) together with the PD. Also plotted are the same histograms shifted by St = —57.8 ns (green online). 



5t 1 =t 1 - ti = -6.7 ns, (17) 
5t 2 = t 2 - t 2 = 44.5 ns, (18) 

with an average of 18.9 ns. 

To compare the shift of St = —57.8 ns with a 
statistical fluctuation of at least the same size, the 
MC generation was repeated one million times and 
in this process 7 608 samples (0.76%) with statis- 
tical fluctuations Si < —57.8 ns were encountered. 
The PD of the first of these MC samples are shown 
in Fig. 5. They feature the time averages 

<i = (5907±34)ns, (model 1) , (19) 
h = (6575 ± 34) ns, (model 2) , (20) 



with the actual statistical fluctuations (using t% = 
5906.7 and t 2 = 6575.0) 

511 = ti-h = -35.5 ns, (21) 

51 2 = t 2 -i 2 = -83.7 ns, (22) 

resulting in an average of —59.6 ns. Figures 6 and 
7 zoom into the leading and trailing parts of the 
departure times of Fig. 4. 

The shift by —57.8 (from red to green in the 
online version of Figs. 6 and 7) exceeds now the 
binsize of 50 ns, while the noise of the histogram 
entries increases. By visual inspection it is hardly 
possible to distinguish a shift of all times from a 
statistical fluctuation. The analysis of the uniform 
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Fig. 7. Model 2: Zoom of the leading (left plots) and trailing (right plots) MC generated departure time histograms (red 
online) together with the PD. Also plotted are the same histograms shifted by St = —57.8 ns (green online). 
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Fig. 10. Enlargement of the left part of Fig. 8. 

distribution in section 2 suggest to display the tails 
of the distribution event by event. Using the em- 
pirical CDF this is done in the next section. 



well suited for a display of the tails of a distribu- 
tion, and based on the binomial distribution we 
will develop their quantitative analysis relying on 
uniformly distributed random variables, similarly 
to the well known goodness of fit [6] . 

The exact CDF of our PD models labeled by k 
are 



F k (t) 



dt'p k (t'), k = l,2, 



(23) 



where p k {t) is given by (10). The empirical CDF 
E k (t) of models k = 1,2 are calculated from the 
data (see, e.g., [6]), which are here MC gener- 

, i£ and 
< tl i+1 holds for 
0, . . . , n + 1, where the definitions tZ° = — oo 
+oo have been added. The empirical 
CDF are then defined by 



ated random times. Let the data be t\, 



r 1 



,t£™ sorted, 



and t£" +1 



5. Cumulative distribution functions 
(CDF) 

For continuous distributions a problem of his- 
tograms is that a binsize needs to be chosen. In 
our case we need a good resolution on the time 
axis, i.e., a small binsize. However, this increases 
the noise of the bin average. As the empirical CDF 
has no such free parameter, it appears more suit- 
able for the analysis at hand. In particular it is 



E k {t) 



n 

0,1, 



for tl* < t < tl' 
, n, n + 1 . 



Shifted by St the empirical CDF are 
S k (t) = E k {t-5t). 



(24) 



(25) 



For the data sets of random times used in the pre- 
vious section F k (t), E k (t) and S k (t) with St = 
—57.8 ns are displayed in Fig. 8. The scale on the 
right ordinates gives the usual probability defini- 
tion for which a CDF is in the range [0, 1], while 
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Fig. 8. Tails of the CDF (F^ exact, empirical, Sk empirical shifted by —57.8 ns, k = 1, 2 model). 
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Fig. 9. Tails of the CDF for departure times with a —59.6ns statistical fluctuation (F^ exact, E^ empirical). 



Q 

O 



the scale on the left ordinates is chosen to count 
the number of included data points. Fig. 8 should 
be compared with Fig. 9, which features a statis- 
tical fluctuation of the average time by —59.6 ns. 
While the — 57.8 ns shift is clearly visible for all four 
cases in Fig. 8, there is no similar signal in three of 
the four cases depicted in Fig. 9. In the figures on 
the right side Fk (t) = 1 is only slowly approached, 
because in our approximation of the PD there are 
still small contributions from the region starting 
between 15 000 ns and 20 000 ns (about 30 events 
for model 2 for which the effect is larger than for 
model 1). These events do not spoil the over- all 
picture, so that the accuracy of a PD in this region 
does not really matter. 

Visually, shift and statistical fluctuation behave 



already differently, though not always. The task 
is now to develop quantitative criteria, which are 
in this section based on the binomial distribution. 
Defining 



Ft = F k {t?) and G\ = l-Fl 



(26) 



the probability for a single data point t to be in 
the range < Fk (t) < Ffc is FJ, and its probability 
to be in the range F l k < F k (t) < 1 is G{. The 
probability that i out of n data points are in the 
< F k (t) < F l k range is 

a - j^m {H)< (G * r * ■ (27) 

If there is a shift by St < we expect that too many 
data points populate the < Fk{t) < FJ, range, so 
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Fig. 11. Average Q % k as functions of i from 10 000 samples. 



that 



(28) 



becomes small. It is easy to see that the Q\ are uni- 
formly distributed random variables, when the hy- 
pothesis is correct that the data points are created 
with the distribution Fk (t) . So the usual interpre- 
tation [6] of a likelihood that the discrepancy be- 
tween the hypothesis and the data is due to chance 
applies to the Q\. 

To present an example we give in Fig. 10 an 
enlargement of the left tail from Fig. 8 by pick- 
ing out the 32 smallest times per model. Upon 
inspection of the numbers one finds for model 1 
tl 1 = 359.985 ns and = F^ 1 ) = 5.7 x 10" 9 . 
This implies Q\ = Ql = 4.3 x 1CT 5 for model 1. 

29.968 ns, which gives F\ = 
3.5 x 10~ 6 and for n — 7612 one finds 
Q\ = 0.026. Combined, these two events alone give 
a probability Q = Q\ Q\ = 1.15 x 10" 6 for the 
likelihood that the discrepancy is due to chance. 

To understand the typical behavior encountered 
for our two distribution functions, we generated 
10 000 MC samples. In Fig. 11 average values for 
Q\ are shown. The almost straight lines at 0.5 are 
the averages obtained when one generates random 
times with the underlying PD. The two other lines 
correspond to the average Q\. obtained when each 
sample undergoes a shift of —57.8 ns. 

At a first glance the average Q\ values for the 
smallest times look far less promising than what we 



For model 2 t^ 1 



mi 1 ) 



Fig. 12. Probabilities for Q\ < p, St = -57.8 ns, from 
10 000 samples. 

may have hoped for from the inspection of our trial 
sample. For model 1 the average is at 0.017 and 
for model 2 even at 0.25. At a second look we shall 
see that this comes mainly from the outliers of the 
non-Gaussian distribution of the Q\, while for a 
majority of cases our trial sample remains typical. 
Besides Q\ for the smallest times, an attractive 
choice for Q\ is at or close to the minimum of the 
average, which we choose to be i = 80 for model 1 
and i — 160 for model 2. Due to the increased 
number of data there are fewer outliers than for 
Q\. Note that for different i values the Q\ are not 
statistically independent, because the calculation 
of Q\ includes all data already used for Q 3 k with 
j < i. 



We define Py.{p) to be the probability for Q\ 



< 



p. Fig. 12 plots Pi (p) and we find Q\ — for more 
than 80% of the samples and Q\ = Q for almost 
23% of the samples. To combine the results from 
both models we notice that 



Q\ Q2 



l-HQ\Qi) 



(29) 



is again a uniformly distributed random variable 
with the interpretation that the discrepancy be- 
tween data and hypothesis is due to chance (Q 1 ^ — > 
Q1Q2 for Q\ or Q 3 2 to zero). We use the notation 
Q l 3 = Q l 3 l and find Q\ = for more than 85% of 
the samples, where (29) is of course performed be- 
fore sorting. 

We define now Qf in = Qf, Q!f in = Qf° and 
Qf in the combination (29) of the two. P™ in (p) are 
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Fig. 13. Probabilities P™ in for Q™ in < p, St = -57.8 ns, 
from 10 000 samples. 
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Fig. 14. Diagonal lines: Probabilities P™ in for Q™ in < p, 
St = ns. Upper curves (F): The same with a statistical 
fluctuation St < —57.8 ns of each sample average. Each 
curve relies on 10 000 samples. 

the corresponding probabilities for Q™ m < p. In 
Fig. 13 the P™ n (p) are plotted versus p using a log 
scale on the abscissa. The Q™ m probabilities tend 
to be small, but never zero. For the combined data 
we find Q^ in < 1(T 7 for about 14% and Qf n < 
1CT 3 for more than 82% of the samples. Now Qf in 
works quite independently from Ql, reducing the 
situations, where their combined estimate Qf ln — 
Qf n Ql [i - HQf n Ql)] gives Qf in > io- 3 to 
approximately 2.6%. Compare Fig. 13. Although 
the application of Eq. (29) is not entirely correct in 
this situation, the related bias is small as it comes 
from the contribution of single data points t^} to 
each of the Qf n . 
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Fig. 15. Probabilities P* for Q\<p, St = ns and (F) the 
same with a statistical fluctuation St < —57.8 ns of each 
sample average. Each curve relies on 10 000 samples. 



For the peace of mind we show in Fig. 14 the 
pmm vames wnen no shift is applied to the data. 
As claimed the Q™ m turn out to be uniformly dis- 
tributed random variables in the [0, 1) range with 
a small bias visible for Q" lm . These are the almost 
straight P™ m (p) = p lines. The upper part of the 
figure shows P™ m {p) from 10 000 samples, which 
are selected so that each exhibits a statistical fluc- 
tuation Si < —57.8 ns (as about 0.76% of the sam- 
ples exhibit such a fluctuation it requires the gener- 
ation of more than 13 million samples) . Clearly, the 
effect on Q™ m is minor and -P™ ln (p) probabilities 
for Q™ m < p are only slightly enhanced compared 
to the straight line. The same analysis shows for Ql 
practically no deviation from the straight (p) — 
p lines as is demonstrated in Fig. 15. These state- 
ments remain true when a bias of the type dis- 
cussed in [4] , which enhances the initial tail of the 
PD, contributes to the shift in the mean value St. 

In summary, there is high probability that a 
study of the initial tails of the departure time dis- 
tributions reveals whether there is a statistically 
relevant shift or not. A similar study could be per- 
formed for the back tail of the departure times. 
One could further expand the analysis of this sec- 
tion to provide an explicit estimate of the shift St, 
similarly as with Eq. (9) for the uniform PD. 
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Fig. 17. Log-likelihood functions (red online) of our typical sample for model 1 and 2. Parabolic fits (green online) are 
indicated for the central regions. 



6. Maximum likelihood method 

For the statistical estimate of the time shift value 
St a maximum likelihood method is the choice of 
Ref. [1]. Its properties are studied in this section 
using 10 000 MC generated samples for each of the 
PD. 

The log likelihoods for our PD (2) are given by 



lni fe (t) = J2]n[p k (f+t)] , k = 1,2 



(30) 



»=i 



where the times f are MC generated with the PD 
Pk (t) . It is instructive to replot Fig. 1 for logarithms 
of the probabilities. This is done in Fig. 16 with 



Pk(t) normalized to Y^iPk(t % ) = 1 where (in ac- 
cordance with our discretization) t 1 is incremented 
in 1 ns steps. These plots exhibit clearly the loca- 
tions of small and zero probabilities, which data 
will (statistically) avoid when they are created with 
the underlying PD pk(t). However, when shifting 
the t % values by t, regions of low probabilities may 
be hit as discussed in the previous sections. 

For our typical sample [7], the resulting ln£fc(t) 
functions are displayed in Fig. 17. For model 1 we 
see for decreasing t, just before — 60 ns, a sharp 
drop ln_Li(t) — > — oo. This comes because the 
smallest time, t 711 + 1 hits a tiny probability, which 
we already encountered in the previous section. 
For model 2 the log- likelihood function lnL 2 (i) 
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Fig. 18. Log-likelihood function (red online) for model 2 
from one of our first ten samples. A parabolic fit (green 
online) is indicated for the central region. 
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Fig. 19. Log-likelihood function (red online) for model 2 
from a sample for which a parabolic fit does not work. 

is smooth down to t = —65 ns, while for positive 
t > 12 ns the largest data point £* n + t hits a 
region of very small probability. 

Fig. 18 shows a model 2 log-likelihood function 
In Z/2 (t) for another of our first ten samples and this 
plot is quite similar to the one for extraction 2 in 
[1]. In the samples one finds a variety of deviations 
from Gaussian shapes, small ones due to the bulk 
structure of the PD and large ones due to a few 
events in the tails. 

Following [1] a step further, parabolic fits are 
made to the central region of each sample. To au- 
tomatize the fitting procedure the central region 
of a model k sample was defined to be the t range 
connected with St k so that ln£fc(i) > lnL™ ax — 10 



holds, where the maximum of the log-likelihood 
function is 

In Lf** = max [In L k (t)] = In L k {St k ) . (31) 

Such parabolic fits are included in Fig. 17 and 18. 

Let us investigate St k estimates for random times 
generated with the PD (10). In the limit of infinite 
statistics 

lim - lni fc (t) = y^Pfe(i) \np k (i + t) (32) 

n— >oo fi * — ' 

i 

holds, where the sums (k = 1,2) corresponds to 
our discretization time steps of 1 ns. They include 
Pk(i) — contributions. For t — they are zero 
due to 

lim (p hxp) = . (33) 

For * ^ 0, e.g. already t — ±1 ns, we can create mis- 
matched contributions Pk(i) mpfc(* + t) with p = 
p(i) finite and q = p(i + t) =0, so that 

lim (p Inq) = — oo (34) 

q— >0 

becomes possible. Hence, in the limit of infinite 
statistics 

lim (6t k ) n = = lim ((^ fe ) 2 \ , (35) 

n— too " n—too \ In 

where (. . .)„ are expectation values for samples of 
n data per model. Deviations of St k from zero re- 
flect statistical fluctuations and bias due to a finite 
statistics. In the following we investigate this for 
n = 7 612. 

Direct estimates of 5t k are obtained by simply 
evaluating the log-likelihood functions (30) for a 
sufficiently large region around t = 0, which is here 
taken to be [—65ns, 65ns]. The positions of the 
maxima of the parabolic fits give also estimates of 
St k , which deviate to some extent from the direct 
estimates. For instance, for our typical sample of 
Fig. 17 we find by direct estimate (t in steps of 1 ns) 

Sti = 5 ns and St% = 6 ns , (36) 

while the parabolic fits give 

<fti = (2.3 ±9.9) ns, 5t 2 = (3.7 ± 10.6) ns, (37) 

where, as in [1], the error bars are standard devi- 
ations obtained from the parabolic fits under the 
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assumption that the center of the likelihood func- 
tion is described by a Gaussian distribution. That 
the Stk values in (36) and (37) are both positive 
is an accident. E.g., for the sample of Fig. 18 one 
finds by direct estimate St-z — —6 ns and 5t 2 — 
(—4.7 ± 11.4) ns from the parabolic fit. 

The analysis of our 10 000 samples allows to im- 
prove on the questionable Gaussian assumption. 
It turns out that the direct estimates of the Stk 
values from lnL™ ax = In Lk(Stk) are more robust 
than the estimates from parabolic fits, which can 
become unstable due to non-Gaussian behavior. 
One rare example is shown in Fig. 19. In that case 
a small t shift, either to the left or to the right, 
does immediately shift one of the extrema t^ 1 or 
ij 71 hito regions of very small likelihood, so that a 
parabolic fit becomes impossible. 

For a finite statistics the extraction of Stk from 
a given log-likelihood function In Lk(St) is a non- 
linear procedure, so that we have to anticipate a 
bias [6] , which means 

(5t k )n + 0. (38) 

The exact maximum position, from which the de- 
viation by a shift t has to be calculated, is at (Stk) n - 
This bias falls off with 1/n, so that it tends to get 
swallowed by the l/y/n behavior of the statistical 
noise. For the direct Stk estimates we find in our 
models 

(<y*i>78ia = (-0-62 ± 0.10) ns, (39) 
(ft 2 ) 7 6i2 = (+0.22 ± 0.11) ns, (40) 

which is in each case smaller than our discretiza- 
tion. The standard deviations were 

ASt 1 = 9.8 ns and ASt 2 = 10.2 ns. (41) 

Combining both models gives the number of the 
abstract 

ASt = 7.1ns. (42) 

To obtain estimates of these numbers from 
parabolic fits one has to eliminate 19 samples for 
which the fits are erratic. Afterwards, averaging 
the standard deviations of the fits gives 

ASh= (9.7±0.9)ns (43) 
ASt 2 = (10.1 ± 1.1) ns, (44) 
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Fig. 20. Empirical peaked CDF of St estimates from 10 6 
samples (red online) in comparison with the peaked Gaus- 
sian CDF of standard deviation 7.1ns (green online). 

where the error of the error is with respect to esti- 
mates on a single sample. Hence, for both models 
combined 

ASt = (7.0±0.7)ns. (45) 

The ASt = 7.5 ns value from the Gaussian fit to 
our typical sample is well consistent. Bias estimates 
stay much smaller than the standard deviation, but 
differ from (39) as the fitting itself imposes a new 
non-linearity. We abstain from giving these num- 
bers, because they reflect also details of the fitting 
procedure like our definition of the central region. 

As the distributions of the Stk are non-Gaussian, 
one may wonder whether Gaussian confidence lim- 
its provide a correct interpretation of the error 
bar (42). Indeed, the empirical peaked CDF [6] 
found for the St estimates via the direct method is 
broader than the Gaussian peaked CDF with stan- 
dard deviation 7.1 ns. See Fig. 20, which relies for 
the empirical peaked CDF on 10 6 random times 
samples. The largest fluctuation in the negative di- 
rection is at —35.5 ns which has thus a probability 
of approximately 10 -6 . This is safely away from 
-57.8 ns. 

Finally in this section, we address correlations 
between the statistical fluctuations of Stk and those 
of Stk- Using i = l,,.,,N — 10 000 samples, we 
denote Stk for sample i by ST k and Stk by (5i™ ax by 
St\. They have expectation values 

(<S?k)78i2 = (St\ - bias) 76 i 2 = . (46) 
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Fig. 21. Scatter plot of (<5l^,<5tj — bias), i = l,...,n, 
N = 10 000 (model 1). 

The bias correction is so small that it can as well 
be ignored. Defining the correlation coefficients by 

C(St kM = %j|M^g) , (47 ) 

their calculation from our N = 10 000 samples 
gives 

(C(5h,5h)) 761 2 = -0A7, (48) 

(C{5t 2 ,6t 2 )) 7612 = -QA3. (49) 

This, at the first glance, non- intuitive anti- 
correlation is for model 1 displayed in Fig. 21. It 
may be explained as follows. When St fluctuates to 
negative values, there is an overpopulation of small 
t values in the sample. Therefore a shift in the 
same direction has a higher probability to produce 
large negative contributions to the log-likelihood 
function (30) than a shift in the opposite direction. 

7. Summary and Conclusions 

We have investigated the problem of identify- 
ing a shift in a broad, non-Gaussian distribution 
of data. For two model probability densities (PD) 
shown in Fig. 1 methods are developed and illus- 
trated, which allow in samples of 7 612 MC gener- 
ated departure times to identify a time shift with 
a precision of about ±7.1 ns while the background 
noise of statistical fluctuations is with ±23.8 ns 
much larger. 



When the boundaries of the distribution are suf- 
ficiently sharp, uniform probabilities from binomi- 
als of the CDF provide excellent indicators. Sub- 
sequently it has been demonstrated that the max- 
imum likelihood method gives a quantitative esti- 
mate of the shift St and its statistical error, while 
bootstrap simulations allow one to avoid Gaussian 
assumptions. Obviously, the analysis of this paper 
can be re-done for other desired sample sizes. 
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